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Abstract 


This  report  comprises  the  final  statement  regarding 
some  initial  exploratory  research  done  by  the  Institute 
for  Acoustical  Research  (IAR)  on  the  Hostile  Artillery 
Location  (HAL)  project.  The  objective  of  this  research 
was  to  re-examine  experimental  data  collected  by  Honeywell, 
Inc.  and  the  Environmental  Research  Institute  of  Michigan 
(ERIM)  during  December  1975.  The  data  collected  comprised 
a series  of  digital  tapes  recording  the  response  of  a 
seismometer  array  to  the  firing  of  various  howitzers  at 
various  ranges. 

IAR's  analysis  of  this  data,  though  not  considered 
technically  complete,  nevertheless  is  considered  very 
successful  at  least  with  respect  to  bearing  estimates 
at  5 km.  The  composite  IAR  estimate  on  three  firings 
is  within  a few  hundredths  degree  of  true  and  the 
standard  deviation  is  1^  degrees.  Range  estimates 
can  be  improved  with  a larger  array  and  a more  so- 
phisticated second  stage.  S 

This  work  was  supported  by  the  Office  of  Naval 
Research  (ONR) , Earth  Physics  Program  (Code  463) . 

The  contract  number  is  N00014-77-C-0446. ' 
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Section  I 


Introduction 

This  report  comprises  the  final  statement  regarding 
some  initial  exploratory  research  done  by  the  Institute 
for  Acoustical  Research  (IAR)  on  the  Hostile  Artillery 
Location  (HAL)  project.  The  objective  of  this  research 
was  to  re-examine  experimental  data  collected  by  Honeywell 
Inc,  and  the  Environmental  Research  Institute  of  Michigan 
(ERIM)  during  December  1975.  The  data  collected  com- 
prised a series  of  digital  tapes  recording  the  response 
of  a seismometer  array  to  the  firing  of  various  howitzers 
at  various  ranges. 

IAR's  analysis  of  this  data,  though  not  considered 
technically  complete,  nevertheless  is  considered  very 
successful  at  least  with  respect  to  bearing  estimates 
at  5 km.  The  composite  IAR  estimate  on  three  firings 
is  within  a few  hundredths  degree  of  true  and  the 
standard  deviation  is  l*s  degrees.  Range  estimates  can 
be  improved  with  a larger  array  and  a more  sophisticated 
second  stage. 

Study  of  the  5 km  events  has  led  IAR  to  devise  a 
processing  scheme  which  is  believed  will  be  successful 
at  11  km  and  possibly  at  17  km.  The  scheme  has  not 
yet  been  implemented. 
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The  processing  employed  on  the  5 km  events  is 
treated  in  great  detail  in  sections  II  and  IV.  Section 
III  supplies  necessary  background  on  the  ERIM  experiment. 
Processing  to  date  on  the  11  km  events  is  presented  in 
section  V.  Section  VI  gives  the  improved  processing  to 
be  implemented  on  the  long  range  events.  Finally  there 
are  a series  of  Appendices  giving  necessary  derivations 
and  technical  background  to  support  the  text. 

This  work  was  supported  by  the  Office  of  Naval 
Research  (ONR) , Earth  Physics  Program  (Code  463) . The 
contract  number  is  N00014-77-C-0446. 


Section  II 


BASIC  PROCESSING 


This  Section  presents  a brief  though  precise  discussion 
of  the  mathematical  procedure  by  which  IAR  has  processed  the 
Honeywell  ERIM  data.  The  mathematics  is  deliberately  simple 
for  this  first  approach  to  the  problem  and  much  has  been 
learned  to  prepare  the  way  for  later  studies.  In  particular, 
the  approach  to  date  has  been  to  regard  the  Haar  transform 
as  a detector  of  local  nonstationary  behavior,  and  this 
seems  valid.  However,  it  has  been  noted  as  a result  of 
these  studies,  that  its  success  stems  from  its  ability  to 
whiten  the  background.  The  Haar  transform  is  a fixed  (in 
time)  linear  operator  and  thus  works  best  on  stationary 
or  quasi  stationary  data.  At  great  range,  however,  the  S/N 
ratio  is  small  enough  that  a detection  system  must  contend 
as  well  with  non  stationary  ambient  noise;  therefore, 
recently  IAR  has  begun  to  focus  on  generalizations  of  the 
Haar  transform  which  adapt  to  non  stationarities  in  the 
background  and  still  whiten  them  relative  to  the  "signal". 

This  broad  class  of  operators  are  the  adaptive  auto- 
regressive operators  and  they  show  great  promise.  As  an 
example,  the  Haar  transform  detects  a low  frequency  transient 
on  7 phones  of  event  39  at  11  km.  This  detection  is  sharpened 
by  these  operators.  This  processing  is  described  in  Section  V. 

The  implemented  IAR  detection  scheme  is  conceptually 
divided  into  two  stages.  Stage  one  obtains  a set  of  arrival 
times  of  transients,  and  stage  two  processes  these  to  deter- 
mine whether  they  comprise  a signal,  and  if  so  where  did 


r i 


-3- 


-4 


it  originate. 

Explicitly  the  stages  are  as  follows: 

First  Stage 

At  a particular  geophone,  a time  series  x(t)  is 
measured.  It  is  then  lowpass  filtered  (anti-alias)  at 
seme  frequency  F and  sampled  appropriately 
Let  the  samples  be  denoted  as: 

X(n)  r Xp(ndt) 

The  discrete  sequence  x(n)  is  organized  into  blocks  of 

M 

length  N (a  power  of  two,  N « 2 ) samples  each. 

Let  x(n):  n = 0,...,  N-l  denote  such  a block,  and  let 
^(k)  denote  its  kth  Haar  coefficient  for  some  sequency 
group  m: 
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The  original  sequence  can  be  retrieved  from  the 
transform  coefficients  as 

C-i  . 

1 2-  2 Xnfo) 

The  Haar  functions  Hm(k,t)  form  an  orthogonal  set  of 
functions : 


/ 0>,t:  HM.t)  At  --  <r„r.  <Fke 


It  is  seen  that  the  Haar  transform  (for  n^O,  k^O) 
effectively  forms  the  difference  between  the  pair  of  two 
mean  values  computed  over  adjacent  subregions  of  the  data 
block.  Each  local  mean  is  computed  over  N/2n  samples. 
Consequently,  each  Haar  coefficient  (differenced  local 
mean)  is  affiliated  with  a data  sub-block  whose  center  is 
located  at  time  t = nAt  = T . The  index  k is 

therefore,  effectively,  a decimated  time  index.  For  m = 0 
the  value  of  the  single  Haar  coefficient  is  the  mean  value 
of  the  time  series. 

The  Haar  transform  therefore  provides  a way  to  examine 
the  change  in  the  mean  value  of  a time  series.  Note,  as 
well,  that  if  the  time  series  is  squared  before  being  trans- 
formed, then  the  coefficients  of  the  transform  will  represent 
changes  in  the  second  moment  of  the  original  series.  These 
operations  could  be  done  by  other  means,  but  an  advantage 
of  the  Haar  transform  is  that  there  exists  a fast  computer 
algorithm  for  the  computation  of  the  coefficients,  based 
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upon  the  factorability  of  the  matrix  Hm(k,t).  For 
sequences  of  length  1024  pts  the  Haar  transform  is 
roughly  35  times  faster  than  the  FFT.  Even  greater 
speed  advantage  is  attained  if  the  digital  computer 
utilized  does  not  have  floating  pt.  hardware,  or  if 
the  sequences  are  longer. 

If  the  input  sequence  is  a sample  of  a stationary 
time  series  which  is  independently  normaly  distributed 
with  known  parameters  N£w,ffA)  then  the  Haar  coefficients 
xm(k)  will  be  independent,  normally  distributed  as 

.l»W  | . 

N(0,  ^ - <J*“  ).  If  the  variance  is  unknown  the  Haar 

coefficients,  normalized  by  the  sample  standard  deviation, 

i s'"*' 

will  be  t-distributed  with  degrees  of  freedom  x1  • ~ - 

We  assume  that  the  ambient  seismological  noise  series 
is  stationary  so  that  the  mean  is  constant.  The  above 
described  t-statistic  will  then  fluctuate  about  zero. 
Significant  departures  from  zero  will  then  be  correlated 
with  local  non  stationarity  in  some  sequency  group. 

(Recall  that  there  is  a time  dependent  t-variable  for 
each  sequency  group) . We  define  such  significant  ex- 
cursions of  the  t variable  as  designating  the  reception 
of  a "transient"  wave  form  within  the  data  sub-block 
affiliated  with  the  corresponding  Haar  coefficient.  If 
the  Haar  coefficient  is  given  by  xm(k*) , then  the  sub- 
block  is  the  one  of  length  2N/2m  centered  at  t = T, 

We  refer  to  this  data  sub-block  as  a "contender  segment". 
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Let  us  now  suppose  that  the  above  procedure  has  been 
carried  out  on  all  geophone  channels.  (Note  that  if  a 
rough  bearing  estimate  is  known,  the  data  on  the  horizontal 
geophone  channels  should  be  rotated  to  this  bearing.  This 
procedure  allows  approximate  separation  of  the  data  into 
components  transverse,  longitudinal,  and  vertical  to  the 
incoming  wave  front,  and  thereby  facilitates  separation 
into  Love  and  Rayleigh  wave  components) . 

We  further  suppose  that  transients  are  detected  on 
several  phones  within  some  time  window.  That  is,  let  T 
denote  a maximum  expected  travel  time  across  the  array 

T*  A/^  „ 


where  A is  the  array  aperture  and  c . is  the  minimum 

min 

expected  seismic  wave  group  speed.  And  let  T^  denote  the 
center  time  of  the  contender  segment  on  the  itil  geophone. 
(These  times  should  be  drawn  from  channels  of  the  same 
type) . 

We  assume  the  times  are  within  the  window;  that  is: 
Tmax  “ max  <V 
Tmin  * min  (Ti> 
then  Tmax-Tmin  < T 


In  this  case,  we  further  analyze  this  set  of  contender 
segments  to  determine  whether  a signal  is  present.  Roughly 
speaking,  a signal  is  present  provided  the  relative  times 
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between  the  transients  allow  the  second  stage  (below) 
to  obtain  spatial  positions  whose  covariance  ellipse  is 
small. 

In  order  to  obtain  a more  precise  estimate  of  the 

relative  times,  cross  correlation  of  the  continuous 

Haar  coefficients  $ (s)  is  carried  out  for  fixed  m 

m 

(m  being  the  sequency  index  whose  t variable  indicated 
a transient) . The  variable  s here  indicates  that  the 
center  of  the  Haar  sub-block  is  varied  continously 
across  the  contender  segment.  This  process  can  be 
visualized  as  Haar  filtration  of  the  contender  segment 
(and  its  adjacent  neighbors).  Briefly  then,  each  con- 
tender segment  within  the  window  is  Haar  filtered  and 
then  cross  correlated  with  a reference  segment. 

Let  the  reference  segment  be  chosen  as  the  center 
phone  of  the  array  (this  is  not  essential) , then  cross 
correlation  of  each  remaining  segment  with  this  reference 
segment  will  give  correlation  diagrams  which  have  a 
maximum  at  some  point.  Let  us  choose  the  relative  time 
between  the  ith  phone  and  the  central  phone  to  be  the 
lag  value  at  which  this  maximum  occurs  on  the  ith  correla- 


tion diagram: 


mn-xR.(r)--  P<  iT}*) 

c 


We  therefore  obtain  relative  arrival  times  between  all 
phones  and  the  central  one.  This  procedure  yields  sharp 


time  estimates  provided  the  correlation  peak  falls  off 

' • 

sharply.  The  precision  of  this  estimate  is  given  ap- 


proximately by  the  formula  (confer  Appendix  A) 

where  S/N  is  the  signal  to  noise  peak  power  ratio  and 
w is  the  bandwidth  of  the  cross  correlated  signals. 

The  bandwidth  can  be  improved  by  carrying  out  the  same 
procedure  on  higher  sequency  groups  on  the  contender 
segment  in  the  manner  of  a "protected  t-test". 

On  the  5 km  ERIM  data,  the  S/N  ratio  is  roughly  16 
and  the  bandwidth  w«  4 Hz.  This  formula  indicates  that 
CTZt  ~ .01,  which  is  consistent  with  the  result 
obtained  from  the  second  stage.  If  the  relative  arrival 
times  are  obtained  by  a procedure  utilizing  jointly  all 
phones  in  the  array  (rather  than  pair  wise  as  here)  this 
precision  can  be  expected  to  improve  as  the  square  root 
of  the  number  of  phones  in  array. 
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Second  Stage 


A set  of  relative  arrival  times  has  been  obtained 


by  the  first  stage,  and  these  times  are  now  processed  by 
the  second  stage  to  determine  if  they  jointly  indicate 
that  their  respective  transients  all  came  from  the  same 
location  in  space.  The  criterion  as  to  whether  this  is, 
in  fact,  the  case  is  given  by  the  area  of  the  x,y  con- 
fidence ellipse.  The  definition  of  the  x,y  confidence 
ellipse  is  that  elliptic  contour  within  which  we  can  be 
95%  certain  that  the  transient  wave  front  originated. 

The  equations  that  follow  explicitly  state  how  the 
x,y  position  is  obtained  from  the  relative  arrival  times, 
and  how  the  confidence  ellipse  is  determined.  Also 
presented  is  the  statistical  means  by  which  the  various 
geophone  channels  are  combined  to  obtain  estimates  which 
are  more  precise.  The  times  T^  below  are  those  obtained 
by  the  cross  correlation  procedure  of  the  first  stage. 
Note  that  the  time  T^  to  the  reference  phone  is  zero. 

The  fundamental  physical  model  is  extremely  simple, 
but  is  successful  on  the  5 km  events.  The  program  is 
not  in  any  way  committed  to  these  equations  and  they  will 
be  made  more  sophisticated  including,  for  example,  dis- 


persion and  inhomogeneity. 

As  depicted  in  Pig.  2.1,  suppose  a source  firing  at  time 
T is  located  at  position  x,y,  and  its  seismic  wave  is  received 
at  times  T^:  i = 1,...,N.  Let  the  ranges  to  these  sensors  be 
denoted  by  r^s  i = 1,...,N.  Then  to  first  order  (assuming 


i |(x.-x)2  +(y . -y) 2 
*- 


These  comprise  n equations  in  4 unknowns  (x,  y,  t,  v)  , and  if 
N = 4 they  can  be  solved  exactly.  However,  the  arrival  times  in 
fact  contain  random  noise  perturbations  which  will  therefore  introduce 
error  into  the  solution.  Consequently,  many  observations  are  made 
and  the  solution  which  maximizes  the  likelihood  function  is  chosen. 


The  model  thus  becomes: 


where  the  Gaussian  random  process  £ has  zero  mean  and  known  covariance 
matrix 


We  choose  x,  y,  T,  v to  minimize  the  form 


Q(x,y,T,v)  = E (T. -h  ) 

i > j 1 1 


where 


denotes  the  inverse  of 


m . . (x)  . 

lj  3xj 


< x< 


1 
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The  normal  equations  for  this  problem  become 


<W 


0 


k*  1 * • • • 1 4 


These  non-linear  equations  are  to  be  solved  for  the  unknown  vector 
estimating  x. 


The  solution  is  obtained  as  the  convergent  limit  of  the  Gauss-Newton 
iteration  scheme 


where : 


> o'*1’  = ‘[(h^c1  H’_1  hT 

X(n)  denotes  the  value  of  x at  the  n iteration 


x(n)  _ t - h (x^) 

T = (TjjTj*  • • • >Tj|) 

;C°)  is  an  initial  guess 

*(0) 

The  iterative  scheme  will  converge  for  virtually  any  choice  o jt 
to  the  same  limit.  In  the  program  implemented  the  matrix  is  taken 

to  be  cr*MC  times  the  identity  matrix. 

KNS 

j HH 
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RMS 


J [Ti  - hi  'i(n,j 


The  estimator  x has  a covariance  matrix  given  by 


- (HT^e‘  H)-1 


The  la  confidence  ellipse  about  ^x  is  derived  from  this  matrix, 


The  above  method  estimates  a velocity  which  should  be  taken  to 
represent  the  mean  horizontal  group  velocity  from  the  source  to  the 
array.  To  a great  extent  the  method  succeeds  on  real  data  because 
this  mean  velocity  is  a very  stable  physical  parameter,  and  its  value 
is  not  assumed  known  by  the  program 


The  equations  governing  the  combination  of  observations  are  as  foil 


ows 


Let  x and  y independently  estimate  a parameter  <x.  Let  A ; and 
denote  their  covariance  matrices.  Then  the  minimum  variance  unbiassed 

A 

composite  linear  estimator  Z is  given  by 


-4.z  • l " A”;  • * « A‘  • X 
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Section  III 


ERIM  Experiment 


This  experiment  was  planned  by  Honeywell,  Inc.  and  the 
Environmental  Research  Institute  of  Michigan  (ERIM) , and 
was  carried  out  during  December  1975.  The  weapon  firings 
were  done  in  the  northwest  corner  of  Anoka  county, 
Minnesota,  and  the  geophone  arrays  were  located  in  Insanti 
and  Sherbourne  counties.  The  relative  locations  are  de- 
picted in  Fig.  3.1. 


The  geology  of  the  test  region  is  given  in  Appendix  c,  but 
it  is  worth  noting  here  that  there  were  a number  of  frozen 
lakes  between  the  firing  site  and  the  geophone  arrays. 
There  was  significant  snow  cover,  but  seme  care  was  taken 
to  ensure  geophone  contact  with  the  ground. 


The  recording  array  consisted  of  nine  3-axis  geophones 
arranged  in  a cross,  as  shown  in  Fig.  3.2.  The  geophone 
spacing  along  each  leg  of  the  cross  was  75  meters.  The 
legs  were  at  right  angles  with  the  north-south  leg  aligned 
with  true  north.  It  is  not  known  how  precisely  this  con- 
figuration was  surveyed. 


The  weapon  was  fired  horizontally  into  a corimgated  steel 
catcher  backed  by  20  feet  of  coarse  sand  between  the 
steel  and  a 100  foot  earth  backstop.  The  azimuth  relative 
to  true  north  was  268.5°  as  shown  in  Fig.  3.1.  All  weapon 
shots  were  a.  charge  7 with  a sand  load. 
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There  were  a total  of  88  firings.  IAR  has  examined  all 
5 km  recorded  weapon  firings,  of  which  three  were  free 
of  significant  hardware  problems  (confer  Appendix  D) . 

These  3 shots  are  fully  discussed  in  section  IV. 

There  remain  77  events  consisting  of  forty  11  km.  events 
and  thirty-seven  17  km.  events.  Appendix  B includes  a 
transcription  of  the  data  logs.  Regarding  these  events, 

IAR  has  taken  the  point  of  view  that  it  is  best  to  con- 
centrate heavily  on  a single  event  at  a given  range  until 
it  is  thoroughly  understood  before  attempting  further 
events  at  that  range.  Moreover  it  is  though  best  to 
understand  thoroughly  the  11  km.  events  than  to  prematurely 
examine  the  17  km.  events. 

There  sure  very  sound  scientific  and  economic  reasons  for 
adopting  this  point  of  view.  The  scientific  reasons  are 
clear,  and  the  economic  reasons  stem  from  the  necessity 
to  obtain  computer  generated  plots  of  the  raw  data. 

This  process  is  very  time  consuming,  and  consequently 
the  IAR  approach  has  been  to  generate  these  plots  until 
an  event  is  found  which  is  free  of  hardware  clipping  and 
electronic  noise.  This  event  is  then  intensively  studied 
as  described  above.  Thereafter  all  other  "hardware  o.k." 
events  at  the  same  range  are  processed. 

IAR  has  examined  the  raw  data  for  the  17  km  events  36,  37, 
38  and  has  selected  event  #38  for  study.  The  progress  to 
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date  on  this  event  has  been  good.  This  matter  is  further 
described  in  section  VI.  The  plotting  of  raw  data  for  the 
17  km.  events  was  chosen  to  begin  at  event  36  because  it 
is  the  first  event  with  a 109  m.m.  howitzer.  This  weapon 
was  chosen  by  virtue  of  the  fact  that  of  the  three  pro- 
cessable  5 km  events  (#8,  9,  10),  two  of  them  were  with 
the  109.  It  should  be  stressed  again  that  "proces sable” 
refers  only  to  freedom  of  the  recording  from  hardware 
problems.  Whether  the  109  m.m.  howitzer  produces  a 
larger  seismic  signal  than  other  weapons,  has  not  yet 
been  examined. 

Appendix  B gives  the  correspondence  between  the  data 
channel  no.  and  the  geophone  - axis. 


Section  IV 


5 Km  Events 

This  section  presents  the  results  of  the  application 
of  the  processing  techniques  described  in  section  II  to 
some  5 Km  events. 

From  Table  B-l  it  can  be  seen  that  a total  of  eleven 
events  were  recorded  at  the  5 Km  range;  of  these  three 
were  TNT  charges  and  two  were  ambient  noise  recordings. 
The  remaining  six  events  correspond  to  the  following 


weapons: 

EVENT  WEAPON 

5 105  mm  Towed  Howitzer 

6 

7 155  mm  Towed  Howitzer 

8 

9 155  mm  self  propelled  Howitzer 

10 

Initial  analysis  of  the  raw  data  for  event  5 showed 
that  a strong  coherent  frequency  (~15  Hz)  was  present 


during  the  recording,  and  since  DC  offset  and  some  clipping 
was  observed  IAR  decided  to  concentrate  efforts  on  a dif- 


ferent event.  Event  6 was  not  analyzed  in  the  belief  that 
it  would  present  the  same  recording  problems  but  recently 
the  record  has  been  plotted  and  it  appears  to  be  a sue- 
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During  the  recording  of  event  7 sane  geophone  lines 
appear  to  have  been  open  since  sente  records  are  characterized 
by  hatfd  clipped  high  frequency  noise.  The  remaining  three 
events  8,  9 and  10  were  successfully  processed  and  a brief 
description  of  the  results  follows. 

Fig.  4.1(a)  represents  the  Haar  spectrum  corresponding 
, to  the  sequency  group  m=5  for  the  N-S  channels  in  event  8. 

A coherent  seismic  signal  is  detected  about  7 seconds  after 
the  shot  is  fired.  Fig.  4.1(b)  is  a plot  of  the  raw  data 
corresponding  to  this  contender  segment  and  Fig.  4.1(c) 
is  the  Haar  filtered  version.  The  filter  used  has  an  impulse 
response  defined  by  the  Haar  function  Hg(o,t). 

Figures  4.2  and  4.3  depict  analogous  processing  for 
the  E-W  and  vertical  channels.  Similarly  event  9 is 
represented  by  Figs.  4.4,  4.5  and  4.6  and  event  10  by 
Figs.  4.7,  4.8  and  4.9. 

The  filtered  contender  segments  (continuous  Haar 
spectrum)  are  then  cross  correlated  with  a reference 
channel  to  determine  the  initial  estimates  of  relative 
arrival  times  (Fig.  4.10).  As  described  in  section  II, 
the  estimates  of  time  differences  can  be  improved  by, 
based  on  the  initial  estimates,  narrowing  and  cross- 
correlating  the  contender  segments  with  higher  sequency 
group  Haar  representations,  but  this  technique  was  not 
followed  for  this  first  report  and  the  results  presented 
are  based  on  initial  time  difference  estimates. 
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The  output  of  the  position  estimator  for  the  three 
events  is  given  in  table  4.1.  For  each  event  two  position 
estimates  are  recorded.  The  second  one  is  a result  of  up- 
dating the  measurement  error  covariance  matrix  with  the 
residual  errors  obtained  from  the  first  estimate.  Although 
simplistic  in  nature  this  method  does  improve  the  bearing 
estimation  for  two  events.  For  the  third  event  the  bearing 
estimate  is  not  in  accordance  with  the  simulation  described 
in  Appendix  E;  in  this  case  refining  the  relative  arrival 
time  estimates  by  crosscorrelating  higher  sequency  Haar 
spectra  would  have  been  in  order. 

Consistently  range  is  underestimated  by  the  position 
estimator.  Initially  this  system  was  designed  for  surface 
waves  and  therefore  the  vector  wave  number  associated  with 
a signal  wavefront  was  assumed  to  be  horizontal.  Given 
the  geological  stratification  of  the  recording  area 
(Appendix  C)  and  the  low  frequency  of  the  wave  train 
detected  it  is  logical  to  think  about  a "refraction 
arrival"  canning  into  the  receiving  array  at  a certain 
critical  angle  (~65°) . The  processing  system  had  no 
capability  of  analyzing  three-dimensional  wave  numbers 
and  the  range  estimated  is  apparently  the  horizontal 
distance  to  the  point  of  radiation  to  the  upper  strata. 
Finally  a composite  estimate  of  position,  following  the 
Gauss-Markov  theorem,  is  given  at  the  bottom  of  table  4.1. 
When  the  three  events  are  considered,  the  bearing  error  is 
negligible.  No  doubt  there  is  an  element  of  luck  here, 
but  still  the  occurrence  speaks  strongly  for  the  technique. 
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Table  4.1 
* 

Position  Estimates  for  5 Km 


True  Range:  R^,  = 5.280  Km 

True  Bearing:  »T  = 301.4° 
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Section  V 


Processing  of  the  11  Km  Data 

This  section  discusses  the  exploratory  processing 
IAR  has  carried  out  on  the  11  km  events.  In  accord  with 
the  philosophy  described  in  section  3,  IAR  has  chosen  a 
single  11  km  event  (#38)  for  this  analysis.  The  goal 
of  this  study  is  to  learn  what  refinements  to  the  5 km 
processing  are  necessary  in  order  to  successfully  approach 
this  regime  of  decreased  signal  to  noise  ratio. 

A comparatively  small  amount  of  time  has  been  spent 
on  the  11  km  events,  and  the  study  is  not  yet  technically 
complete.  However,  some  important  conclusions  have  been 
drawn.  Basically  these  conclusions  are  that: 

(1)  The  autonomous  Haar  operator  must  be  generalized 
to  a time  dependent  operator  capable  of  tracking 
the  ambient  noise  statistics. 

(2)  The  pairwise  cross  correlation  procedure  must  be 
replaced  by  truly  joint  array  processing . 

The  processing  extensions  demanded  by  these  conclusions 
are  described  briefly  in  this  section,  and  in  greater 
detail  in  section  6. 

The  processing  of  section  2 relies  upon  the  Haar 
transform  to  detect  transients  through  their  alteration 
of  the  local  mean,  and  at  that  range  (5  km)  the  procedure 
appears  successful.  Application  of  the  Haar  transform  to 
the  11  km  data  however,  indicates  that  this  feature  is  not 
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sensitive  enough  to  serve  as  a reliable  detector.  Fig. 

5.1  shows  the  application  of  the  sequency  4 Haar  transform 
to  the  raw  data  of  event  #38.  It  is  clear  that  the  number 
of  false  alarms  is  too  great.  (It  should  be  noted  however, 
that  the  11  Jon  air  wave  remains  processable  by  the  Haar 
method) . 

A signal,  however,  is  seen  to  appear  in  the  region 
bounded  approximately  by  the  vertical  lines.  This  region 
occurs  at  about  13  secs  after  firing  and  corresponds  to 
a group  speed  of  about  850  m/sec.  The  zero  crossing  rate 
indicates  an  instantaneous  frequency  of  about  4 Hz,  and 
according  to  Appendix  C this  wave  could  therefore  contain 
Love  and  Rayleigh  components.  The  vertical  channels 
exhibit  significant  correlation  (esp.  channels  6 and  27) 
which  tends  to  strengthen  the  possibility  of  Fayleigh  wave 
(denoted  "R-wave)  activity.  The  strongest  signals  at  5*  km 
are  thought  to  be  Love  waves  (denoted  "G-wave") , and  this 
11  km  possible  R wave  therefore,  indicates  a further  in- 
vestigation of  the  5 km  vertical  channels  may  sharpen  the 
results  the^e.  The  raw  data  corresponding  to  this  Haar 
plot  is  presented  in  Fig.  5.2.  The  region  within  the 
vertical  lines  of  Fig.  5.1  corresponds  to  about  18.5  sec 
on  this  Fig.  5.2.  Again  it  is  seen  that  packets  are 
arriving  in  this  interval  (clearest  in  channels  13,  14, 

15). 

In  the  event  these  disturbances  are  howitzer  induced 
R-waves  the  processing  gain  would  be  improved  by  rotation 
of  the  data  to  the  approximate  bearing  followed  by  joint  • 
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processing  of  the  longitudinal  (with  respect  to  the  wave 
front)  and  vertical  channels.  The  transverse  channel 
would  contain  G-wave  activity.  IAR  has  not  carried  out 
this  operation  because  final  justification  of  the  R-wave 
hypothesis  would  require  that  more  11  tan  events  be 
surveyed,  in  order  to  eliminate  the  possibility  of 
seismic  noise  generated  locally  to  the  array.  It  should 
be  noted  however,  that  the  IAR  propagation  model  is  most 
accurate  on  Rayleigh  waves.  For  this  reason  the  issue 
deserves  further  exploration. 

This  observed  R-wave  possibility  has  led  IAR  to 
examine  the  possible  existence  of  higher  frequencies  in 
the  data.  The  rationale  here  is  that,  depending  on  the 
geology,  R-waves  may  travel  shorter  paths  than  G-waves 
of  the  same  frequency  and  consequently  may  be  attenuated 
less.  If  this  is  so,  the  increased  bandwidth  (expected 
at  higher  frequencies)  would  allow  more  precise  time 
estimates.  In  order  to  examine  this  issue,  a series  of 
power  spectra  were  obtained  using  the  Burg  maximum  entropy 


(auto  regressive)  method.  Spectra  were  computed  by  this 
method  on  successive  512  point  blocks  of  the  raw  data. 
The  auto  regressive  model  employed  64  coefficients  which 
is  sufficient  to  give  good  spectral  resolution  for  fre- 
quencies above  2 Hz. 

The  results  of  this  processing  are  illustrated  in 
Figs.  5.3  through  5.29.  These  are  graphic  reproductions 
of  computer  generated  printer  plots  which  are  too  bulky 

■ - 
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to  include  here.  The  originals  are  on  file  at  IAR.  The 
dashed  curves  denote  ambient  noise  spectra  before  the 
howitzer  firing.  (In  each  case  the  block  chosen  is  #3). 

The  solid  curves  denote  spectra  after  the  howitzer  firing 
corresponding  to  a portion  of  the  signal  region  in  Fig.  5.1. 
(In  each  case  the  block  chosen  is  #13  beginning  about  18 
secs,  after  firing).  Figures  5.3-5.11  are  vertical  channels; 
5.12  - 5.20  are  E-W  channels;  and  5.21  - 5.29  are  N-S 
channels. 

A number  of  observations  can  be  made  from  these  spectra: 

(1)  Signal  is  appearing  on  all  phones  and  all  channels 
in  the  frequency  regions  3-5  Hz  and  20-25  Hz. 

(2)  The  ambient  is  also  strong  in  the  region  3-5  Hz. 

(3)  The  ambient  is  weak  in  the  region  20-25  Hz. 
therefore  the  SNR  ratio  is  good  in  this  region. 

(4)  The  vertical  channels  appear  especially  strong  in 
the  high  frequency  region  although  the  peak  fluc- 
tuates in  frequency  somewhat  more. 

These  observations  lead  to  the  following  conclusions: 

(1)  Rayleigh  waves  are  probably  being  received  at 
approximately  20-25  Hz. 

(2)  Simple  processing  methods  may  be  effective  on 
these  waves  due  to  the  good  SNR. 

(3)  The  region  3-5  Hz  shows  premise  but  will  require 
sophisticated  processing  to  extract  the  signal 
from  the  noise.  These  methods  may  be  required 
to  utilize  jointly  (rather  than  pairwise)  all 
array  data;  may  require  adaptive  whitening;  and 
possibly  joint  processing  with  the  Rayleigh  wave 
at  20-25  Hz. 
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The  processing  visualized  as  appropriate  for  the 
low  frequency  band  is  given  in  section  6.  IAR  has  begun 
to  implement  this  processing. 

Before  closing  this  section,  it  should  be  remarked 
that  the  spectra  of  Figs.  5.3  through  5.29  are  almost 
certainly  not  as  sharp  as  could  be  obtained.  The  spectra 
are  computed  frcro  certain  arbitrary  choices  of  block  size 
and  bin  width.  The  processing  of  section  6 is  largely  in- 
dependent of  such  choices.  Moreover  at  high  frequencies 
then  remains  the  possibly  of  poor  ground-geophone  coupling. 

Finally  we  present  a Fig.  5.30  which  is  pertinent  to 
section  6.  Classical  theory  of  detection  in  colored  noise 
entails  whitening  the  background  spectrum  before  performing 
any  correlation  analysis.  Since  the  ambient  spectrum  is  not 
known  a priori,  and  in  any  event  varies  somewhat/  the  most- 
appropriate  prewhitening  technique  is  the  adaptive  auto- 
regressive technique.  (This  technique  also  provides  the 
foundation  of  the  above  mentioned  spectral  line  tracking 
to  be  presented  in  section  6).  For  future  reference. 

Figs.  5.30a  and  5.30b  illustrate  the  effect  of  this  adap- 
tive whitening  operator. 

A nine  coefficient  operator  is  employed  here  and  is 
not  expected  to  be  able  to  whiten  the  low  frequency  lines. 
However,  it  can  be  seen  that  the  filter  does  render  this 
data  more  nearly  white.  Theory  indicates  that  any  pure 
line  will  not  readily  be  whitened  by  this  filter.  Also 
the  adaptation  time  of  the  filter  has  been  chose  to 
correspond  to  the  time  scale  on  which  the  ambient  spectrum 


is  changing,  and  not  to  that  on  which  tha  signal  transient 
makes  its  appearance.  Therefore,  it  is  to  be  expected  that 
this  filter  will  successfully  sharpen  the  results  of  the 
correlation  analysis  on  any  data  containing  packets  of  small 
spectral  width. 
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Fig.  5.2 
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Section  VI 

Processing  Extensions 


This  section  presents  a brief  description  of  the 
processing  IAR  intends  to  employ  on  the  11  and  17  km. 
events.  This  processing  is  a natural  generalization 
of  the  Haar  method  and  will  be  equally  effective  on 
the  5 km.  events.  Moreover  it  is  expected  to  succeed 

on  nontransient  signals  as  well. 

As  before,  the  processing  presented  below  is  nothing 
more  than  a good  point  of  departure  for  the  analysis  at 
the  longer  ranges.  It  has  been  devised  from  experience 
gained  in  the  5 km.  analysis,  and  it  is  expected  that 
changes,  if  needed,  will  be  implemented  during  the 
course  of  study.  In  particular,  the  prdcessing  presented 
here  will  confine  itself  purely  to  the  so-called  "first 
stage"  in  which  arrival  times  are  derived  from  the  raw 
data. 

Recall  that  the  second  stage  position  estimation 
equations  are  to  be  improved  as  well,  although  this 
modification  will  not  be  presented  here  as  it  seems 
secondary  to  the  issue  of  signal  detection. 

Due  to  the  dispersive  and  inhomogeneous  character 
of  the  seismic  medium,  the  temporal  duration  of  a 
transient  wideband  signal  increases  with  range.  The 
phase  characteristics  at  a given  frequency  are  unsteady. 
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There  is,  as  well,  geometric  and  absorptive  attenuation 
(of  which  the  later  is  markedly  frequency  dependent) , 
and  consequently  the  amplitude  decreases  with  range. 

It  has  been  observed  in  the  ERIM  data  that  the 
local  seismic  noise  characteristics  vary  significantly 
from  geophone  to  geophone,  and  at  the  same  phone  there 
is  considerable  spectral  variation  between  even  the  two 
horizontal  channels.  This  ambient  noise  is  not  stationary, 
not  white,  and  not  known  a priori.  However,  a given 
vertical  geological  layering  will  transmit  certain  fre- 
quency bands  better  than  others  and  the  energy  both  in 
the  ambient  and  the  received  signal  is  concentrated  in 
these  bands  though  with  significant  fluctuation. 

The  time  domain  processing  below  has  the  following 
general  structures 

(1)  The  noise  is  independently  (on  each  phone) 
whitened  by  an  adaptive  filter. 

(2)  The  whitened  outputs  are  linearly  combined  under 
the  assumption  that  a plane  wave  of  unsteady 
frequency,  phase,  and  amplitude  is  traversing 
the  array. 

(3)  The  instantaneous  line  spectrum  of  the  joint 
signal  is  examined. 

■ (4)  If  the  "Q"  of  the  line  spectrum  exceeds  a 

threshold  the  local  neighborhood  is  declared  to 
be  a "contender  segment". 
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(5)  All  joint  cross-correlations  of  them 

where  P is  the  # of  array  phones)  are  computed. 

(6)  These  cross  correlations  are  passed  to  the 

p 

second  stage  which  combines  the  j estimates  by 
the  Gauss-Markov  theorem. 

Let  xp(n)  = xp(nAt)  denote  the  raw  data  recording  on 
the  pth  geophone  p = 1,...,P  where  P is  the  # of  Dhones 
in  the  array,  and  let  yp(n)  denote  the  adaptively  whitened 
series.  Then  we  have 

_L 

yp  (n)=  Xp(n)  - 2 CLpgOl)  Xp(n-l ) 

and 

cLf,t  (n)  = cip<(n-i)  V P(n) 

where  Lf/^  are  chosen  to  model  the  seismic  noise  (confer 
Appendix  G) . 

Now  let  -p 

(a)  = /£  /P(n)  cos  (k, /P  + Vy) 

where  (x  , y ) are  the  coordinates  of  the  pth  phone,  and 
P P 

*?  i = 

ki*Sr  t'rn£- 

A 

where  Xj  9 are  wavelength  and  bearing  of  an  (irregular)  plane 
wave.  In  the  data  studies,  we  will  fix  X and  allow  6 to 
vary . although , in  fact,  9 is  approximately  known.  The 
parameter  X can  be  determined  in  the  second  stage,  however, 
this  is  considered  a secondary  issue. 
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We  now  model  the  joint  series  as  an  autoregressive 
series:  Let  ^ (n)  denote  the  predicted  valqe  of  Z(n), 


then: 


jl 


$(Tl)  - Zrn-rn) 


where 


O' 


>r> 


^ ^rr, (r. ) i 2 z(r‘ g (fi )j 


The  parameters  M, r ^ are  chosen  to  model  the  signal 
process. 

The  instantaneous  autoregressive  power  spectrum  at 
time  t = n£t  and  frequency  f is  given  by 


6?  (f  in)  - 


2 At  <r2 


lHCfjn)l 


where  C*  is  the  noise  power  in  the  joint  output  (estimable 
by  recursive  average) . 

denotes  the  complex  function: 


-M 

H(f  ,n)  - d ~ Z e*piZIT  ^k] 

4 J 


(confer  appendix  G) . 

Rather  than  compute  the  spectrum,  it  is  faster  and 
for  our  purposes  more  precise  to  examine  the  roots  of 


P(r,i i)»  I ~ t <*J»>  * 


rr\  1 « 


j 
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The  relation  between  P and  H is  clearly: 


However,  it  seems  sufficient  to  define  the  "Q" 

(not  precisely  the  Q of  RLC  circuitry)  as  follows: 

Let  f^,  f2  be  chosen  such  that  for  all  f £ 
the  condition  |P(f,n)j<  £ . The  parameter  £ is  to 

be  chosen  according  to  a statistical  criterion  described 
below.  And  let  |p(f,n)j  be  minimized  (within  this  in- 
terval) at  fQ.  Then  define 


Qw* ,n) " 


(This  corresponds  to  a power  weighted  Q) . 

If  Q (f Q,n)  > T for  a number  of  time  steps  nAt  then 
the  temporal  region  is  declared  to  be  a "contender  segment". 
The  threshold  T is,  like  £ , to  be  chosen  as  a result  of  data 
studies.  The  two  thresholds  together  govern  the  probabilities 
of  detection  and  false  alarm.  Their  proper  choice  is  based 
upon  the  statistics  of  Q.  If  this  technique  is  successful 
then  these  statistics  will  be  explored  and  rational  choice 
of  £,  T will  be  decided.  In  an  operational  environment 
they  will  be  adaptively  selected  based  on  these  studies. 

The  remaining  parameters  above,  namely,  X.  * ^ j 
are  not  critical  (cf  Appendix  G) , and  are  intuitively 
clear  being  respectively  the  adaptation  tisies  end  filter 
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Given  the  contender  segments  the  procedure  described 
in  steps  5t  16  above  is  carried  out.  These  steps  depart 
from  the  5 km  procedure  only  in  that  many  more  cross 
correlations  are  carried  out.  This  procedure  is  known 
to  be  optimal  (rf.  Hahn,  W.R. , JASA,  Vol.  58,  No.  1, 

July  1975). 
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Section  2 discusses  the  use  of  the  Hear  transform 
to  isolate  "contender  segments"  on  each  geophone.  Each 
such  segment  contains  a transient  suspected  to  be  a signal 
from  an  event  such  as  a howitzer  firing.  In  order  to 
decide  whether  the  transient  is  a meaningful  signal,  the 
second  stage  of  the  processing  utilizes  the  relative 
arrival  times  to  compute  a position  estimate  with  an 
associated  error  ellipse. 

The  relative  arrival  times  are  obtained  from  pair- 
wise cross  correlation  of  the  filtered,  measured  time 
series  at  each  geophone  with  a common  reference  phone 
(chosen  arbitrarily) . An  equation  giving  the  precision 
with  which  the  arrival  times  can  be  estimated  is  presented 
in  section  2.  This  appendix  presents  a derivation  of  that 
formula. 

Let  yQ(t)  be  the  filtered  time  series  measured  at  the 
reference  geophone,  and  let  y^(t)  denote  the  filtered  series 
on  the  second  phone.  All  that  is  required  is  the  arrival 
time  difference  so  we  can  write  without  loss  of  generality 
that: 

yQ(t)  * S(t)  + nQ(t) 

yx(t)  = S(t-T)  + nx(t) 

where  T (possibly  negative)  is  the  relative  time  by  which 
the  received  signal  at  geophone  1 lags  that  at  geophone  0. 

In  the  above,  S(t)  denotes  the  unknown  filtered  signal 
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! 


* 

i 
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waveform.  All  that  is  assumed  about  this  "signal"  is 
that  it  lies  within  the  contender  segment.  More  pre- 
cisely , within  the  contender  segment  is  assumed  to 
exist  a transient  whose  filtered  forms  on  each  geophone 
have  similar  components.  (Some  disimilar  residue  in  the 
signal  may  be  combined  with  the  noise  so  long  as  it  does 
not  greatly  affect  our  assumptions  regarding  the  noise) . 
Initially  we  assume  the  noise  processes  nQ(t)  and  n^t) 

to  be  stationary,  uncorrelated,  zero  mean  Gaussian  pro- 

2 2 

cesses  with  variances (j*  andO"  respectively. 

no  n-L 

We  choose  that  delay  X which  maximizes  the  a - 
posteriori  probability  density  of  X given  the  measured 
series  yQ,  y1? 

Choose  X such  that  p(X  f yQ,y^)  is  maximized. 

We  have  the  following 

piX/y^y^  = p(yjJXfy0)  p{yQ,*t)/p(y0,y1) 

We  assume  p(yQ,T)  * p(yQ)  pCO 
Therefore 

pCTty^y^  = kx  pm  pCyJr,  y0>  1 

where  k,  is  an  irrelevant  constant.  We  further  asstime  we 
are  maximally  ignorant  of  T,  and  take  it  to  be  uniformly 
distributed. 

Therefore  we  obtain 

pCT^o^)  - k2  p(y1|T,y0) 


- 
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as  the  quantity  to  be  maximized  with  respect  to  X. 

We  note  that  we  can  write 

yl(t)  = yo(t_T)  + n(t) 

where  n(t)  * nQ(t-T)  + n^t)  is  a stationary  zero  mean 
Gaussian  process  with  variance 


(■( 


Let  N denote  the  one  sided  power  spectral  density  of  the 
noise  process  nQ  and  and  let  W denote  the  positive 
bandwidth.  We  have  therefore 


. 2 - Z N W * 2 H (|trc)  • N/4t 


using  the  Nyquist  theorem. 

The  series  Z(t)  ~ y((t}-  ycit  -'£) 
is  zero  mean  Gauss ian  with  variance  CT^Z. 


Consequently  we  can  write 


p(y.  I y0jr)  = exp  [-^/  yuft-T)]2^^] 


Expanding  the  integral  above,  and  accumulating 
constants  of  no  interest  we  obtain 


P(rrl  y„,Y,)  - kA  «P[£  j V,(t>  y,(t-rjeteJ 

'i 


In  this  formula  it  has  been  assumed  that  the  contender 
segment  length  Tj-T.,  has  been  chosen  long  enough  that 


is  essentially  independent  of  X.  The  IAR  procedure  com- 
putes the  exponent  above. 

The  estimated  relative  time  T is  taken  to  be  the  value 
of  X which  maximizes  the  quantity  above  (for  which  it  is 
sufficient  to  maximize  the  exponent) . 

Let  q(T)  denote  the  exponent 


Then 


The  second  and  third  integrals  above  are  assumed 
to  be  approximately  equal,  and  the  fourth  integral  is 
assumed  to  be  approximately  zero  since  the  series  are 


In  the  absence  of  noise  h(1r)  = 0 and  the  maximum 


of  q{T)  occurs  at  T = "X* . We  are  assuming  that  h (X)  is 
small  compared  to  g etr)  near 'IT  = 'tT* , therefore  we  can  assess 
the  extent  by  which  the  maximum  is  shifted  as  follows: 


At  the  maximum 


r n 
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and  <^'(r*)-C> 


therefore 


y'cri  * $"(?*)(£•  f) 


It  follows  that 


T-V* 


Let  A denote  T - 'f* , then 


<^2=  rs"<'t*0  ^ch’ix)2] 


where  denotes  the  expected  value  operator.  We  now 

r-  r,- 


consider  the  term  Ex[hcZ)  J 


Exltirt-f]  - ~L  j f s V«-r)  s’(r-r)  Er[ not)  n(y)l 


where  the  limits  have  been  taken  to  be  infinite  by  the 
earlier  remark 

But  :*Cr>(u.>n(iO]  * n cf(u.-\ /) 

where  is  the  Dirac  delta  function. 


Therefore 


"j&Ch'tf/j.  h 


ffj  [sVuijVa. 


(energy) 


(rms  bandwidth) 


The  quantity  g"  Ct*)  is  computed  as  follows 


Again  T,-T.  is  large  enough  to  include  all  of  Sf  therefore 
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Appendix  B 
Table  B— 1 (cont) 

Transcription  of  data  recording  logs 


after  event  41 
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Appendix  B 
Table  B-l  (cont) 


Transcription  of  data  recording  logs 
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Standa  rd 
Position 

5 km  Site 
Recording 
Channel 
(N,  E,  V) 

11  km  Site 
Recording 
Channel 
(N,  E,  V) 

17  km  Site 
Recording 

Channel 
(N.  E,  V) 

1 

13,  14,  15 

1.  2,  3 

1.  2,  3 

2 

16,  17,  18 

4.  5,  6 

4,  5,  6 

3 

19,  20,  21 

7.  8,  9 

19,  20,  21 

4 

22,  23,  24 

10,  11,  12 

22,  23.  24 

5 

7.  8,  9 

13,  14,  15 

13,  14.  15 

6 

10,  11.  12 

16,  17,  18 

16,  17,  18 

7 

1.  2,  3 

19,  20,  21 

7,  8,  9 

8 

4,  5,  6 

22,  23,  24 

10,  11,  12 

9 

25,  26,  27 

25,  26,  27 

25,  26,  27 
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Appendix  C 


MODEL  OF  SEISMIC  WAVE  PROPAGATION  AT  THE  HONEYWELL  TEST  SITE 

From  seismic  ray  theory  and  wave  theory  applied  to  a model  of  seismic 
layering  at  the  Honeywell  Ordnance  Proving  Ground  (HOPG) , we  can  construct 
a graph  of  frequency  vs.  arrival  time  and  phase  velocity  (apparent  wave 
velocity)  vs.  arrival  time  for  the  various  seismic  phases  that  cross  the 
geophone  array:  P waves,  S waves.  Love  waves,  Rayleigh  waves  and  the  sound 
(air)  wave.  This  graph  Is  a guide  to  Identification  of  phases  on  seismograms 
from  the  HOPG  test  program. 

We  first  define  the  model.  The  test  program  at  the  HOPG  Included  shots 
fired  to  a moveable  array  set  up  alternately  at  three  ranges:  5.28  km, 

11.04  km,  and  17.12  km.  Short  refraction  profiles  at  several  points  served 
to  establish  the  nature  of  shallow  layering  in  the  area  (see  Larson  et  al.,  1976) 
From  information  in  that  document,  a representative  geologic  and  seismic 
velocity  section  for  the  area  is: 


- • 

Layer  Thickness 

Compressional 

Velocity 

Shear 

Velocity 

Density 

Post-glacial  outwash 

9.144  m (30  ft) 

762  m/s  (2500  ft/s) 

440  m/s 

1.9  g/ari 

Ice-compacted  glacial  till 

82.296  (270) 

1829  (6000) 

1056 

2.3 

Cambrian  sandstone 

3353  (11000) 

1935 

2.4 

(See  discussion  In  IAR  Interim  Report,  January  1978). 

In  addition  to  simplifying  the  original  refraction  data,  this  table  Involves 
several  assumptions:  that  shear  velocity  is  related  to  the  compresslonal 
velocity  observed  in  the  refraction  profiles  by  Poisson's  ratio  = 0.25  or 
congressional  velocity  = fT  x shear  velocity;  representative  densities  for 
similar  materials  are  taken  from  Clark  (1966)  in  absence  of  values  measured 
at  HOPG;  that  layers  are  homogeneous  and  have  plane,  parallel  upper  and 
lower  boundaries. 
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For  this  model  Me  have  calculated  the  arrival  times  for  congressional  and 
shear  body  waves  (P  and  S waves)  refracted  through  the  third  layer 
(at  3353  m/s  or  1935  m/s).  This  represents  the  least  time  path  for  P and 
S waves  at  all  three  ranges,  assuming  that  there  is  no  material  of  higher 
velocity  at  depth  beneath  the  sandstone.  Actually,  the  presence  of  granite 
at  a depth  of  about  1 km  may  cause  the  body  waves  to  arrive  substantially 
earlier,  especially  at  11  and  17  km,  than  predicted  by  the  model  of  shallow 
layering  that  was  used  for  calculation. 

Me  also  have  calculated  the  group  arrival  time  vs.  frequency  and  the 
associated  phase  velocity  for  Love  and  Rayleigh  waves,  assuming  shallow 
layering  as  In  the  above  table.  This  calculation  was  done  by  the  method 
of  Haskell  (1953).  Results  are  shown  in  the  figure  below.  These  results 
are  affected  very  little  by  the  granite  or  deeper  layering  because  Love 
and  Rayleigh  waves  are  guided  waves,  with  motion  confined  mainly  to  the 
surface  wave  guide.  In  the  lower  part  of  the  figure,  group  arrival  time 
is  shown  as  a function  of  frequency.  These  curves  represent  the  frequency 
of  the  wave  motion  that  is  seen  on  the  recording  at  various  times  measured 
from  the  time  of  the  shot.  There  are  separate  curves  for  Love  and  Rayleigh 
waves  at  each  of  the  three  ranges  of  interest.  The  upper  part  of  the 
diagram  shows  the  phase  or  wave  velocity  associated  with  each  frequency 
and  arrival  time  represented  in  the  lower  part.  This  is  the  velocity  of 
individual  waves  crossing  the  array  as  measured  by  correlation  of  peaks 
or  wave  fronts  from  one  geophone  to  the  next.  Finally,  the  figure  shows 
the  arrival  time  of  a direct  air  wave  travelling  at  343  m/s. 

The  air  wave  and  the  P and  S waves  are  represented  as  undispersed, 
l.e.,  all  frequencies  arriving  at  the  same  time.  On  the  other  hand,  the 
Love  and  Rayleigh  waves  are  dispersed  due  to  the  velocity  layering  of  the 
wave  guide:  lower  frequencies  or  longer  waves  generally  arrive  earlier 
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than  higher  frequencies  or  shorter  waves.  There  are  also  some  bands  of 
Inverse  dispersion  where  lower  frequencies  arrive  after  higher  frequencies. 
These  are  due  to  the  sharp  velocity  discontinuities  between  layers.  If 
the  boundaries  are  gradational,  or  If  there  are  velocity  gradients  rather 
than  discontinuities,  then  the  dispreslon  curves  are  stralghter  than  for 
this  model,  and  may  lack  Inverse  portions. 

In  addition  to  the  simplified  approximations  of  layering,  the  model 
may  also  depart  from  actual  conditions  In  the  sense  that  there  may  be 
substantial  local  variations  In  the  thickness  of  compacted  till  and 
unconsolidated  glacial  outwash  which  are  not  represented  In  the  model. 

Due  to  these  possible  complications,  the  arrival  times  and  frequencies 
represented  In  the  figure  are  only  approximate  indications  of  the  relation- 
ships which  may  exist  in  seismograms. 

Also,  the  contents  of  seismograms  depend  on  excitation  and  attenuation 
characteristics  not  represented  In  the  figure.  If  a given  frequency  is 
generated  at  the  source  and  propagates  without  strong  attenuation  then 
It  will  appear  on  the  seismogram  at  a time  shown  approximately  by  the 
group  arrival  curves  of  the  figures.  Lower  frequencies  propagate  farther 
with  less  attenuation  than  high  frequencies  In  a given  seismic  medium; 
and  surface  (Love  and  Rayleigh)  waves  spread  two-dlmenslonally  rather 
than  three-dimensional ly  as  do  body  waves.  Also,  recoil  sources  are 
relatively  rich  In  low  frequencies  by  comparison  with  tamped  explosions 
of  similar  energy.  These  factors  suggest  that  low  frequency  Love  and 
Rayleigh  waves  should  be  particularly  useful  for  detection  and  Identification 
of  recoil  sources  at  relatively  large  ranges. 
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Figure  Caption.  Group  arrival  time  vs.  frequency  and  corresponding  phase 
velocity  for  seismic  waves  propagating  on  the  layered  seismic  velocity 
model  given  In  text.  The  lower  portion  of  the  figure  gives  group  arrival 


time  vs.  frequency  for  compresslonal  waves  (P),  shear  waves  (S),  direct 
air  waves  (A),  Love  waves  (L)  and  Rayleigh  waves  (R).  ■ 344  km/sec  Is 

used.  P,  S and  A are  undispersed;  L and  R are  dispersed.  Families  of 

curves  are  given  for  each  of  three  distances  5.28  km,  11.05  km  and  17.12  km. 

* 

The  upper  portion  of  the  figure  gives  the  phase  velocity  associated  with 
the  particular  group  arrival  time  and  frequency  given  In  the  lower  part. 
Phase  velocity  Is  the  velocity  with  which  Individual  waves  cross  the 
local  array. 
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Appendix  D 
ERIM  Data 


(0)  The  ERIM  data  consists  of  a set  of  six  digital  tapes 
containing  seismometer  data  recorded  by  an  array  of  nine 
3-component  seismometers.  The  array  was  responding  to  the 
recoil  of  howitzer  firings  at  ranges  of  5,  11,  and  17  kilo- 
meters. A diagram  of  the  array  is  given  in  Fig.  3.2 
for  a typical  (5  km)  event.  There  are  a total  of  88 
events . 

(U)  The  data  were  low-pass  filtered  at  80  Hz.,  sampled 
(non-simultaneously)  at  300  samples/sec.,  and  then 
digitized  at  14  bits/sample. 

(U)  This  experiment  was  planned  by  Honeywell,  Inc.  and 
ERIM,  and  presumably  the  hardware  was  configured  adequately 
for  the  processing  they  envisioned.  However,  from  the  point 
of  view  IAR  must  employ,  the  following  are  deficiencies  in 
the  data: 

1.  The  array  aperture  is  too  small  to  obtain  acceptable 
range  accuracy. 

2.  The  sampling  technique  was  deficient  with  respect 
to  rms  time  error  across  the  array. 

3.  The  gain  settings  at  5 km  were  18  dB  higher  than 
those  at  11  km.  and  17  km.  In  fact,  the  5 km  gain  settings 
are  so  high  that  they  induce  frequent  non  linear  clipping 
with  consequent  power  supply  drain.  It  is  thought  that 
this  must  lead  to  severe  transient  cross  talk  between 
seismometer  traces. 


4.  The  seismometer  channels  (NS,  EW,  and  vertical 
axes)  themselves  exhibit  low  separation  (perhaps  a result 
of  the  mounting)  so  that  the  large  transverse  Love  wave 
component  appears  as  am  artefact  in  the  longitudinal  and 
vertical  channels  and  masks  those  signals. 

5.  In  order  to  evaluate  our  processing,  some  estimate 
of  geophone  location  uncertainties  is  needed  whereas  none 
is  supplied. 

6.  The  polarity  is  reversed  on  the  E-W  channels  of 
certain  geophones. 


Appendix  E 
Error  Propagation 

In  order  to  study  how  measurement  errors  propagate 
through  the  position  estimator  system  for  the  ERIM  array 
a Monte  Carlo  simulation  was  made.  The  geometry  considered 
is  depicted  in  Pig.  El.  The  source  of  the  perturbation 
was  located  at  ranges  of  5 and  11  km  and  a bearing  of  -1 
radian  with  respect  to  a cartesian  system  defined  by  the 
array  axis.  The  seismic  wave  paths  were  assumed  to  be 
straight  lines  and  the  propagation  medium  was  homogeneous. 
Exact  travel  times  T^e  to  each  sensor  i were  computed  and 
were  perturbed  with  a random  noise  component  T^  obtained 
from  a zero  mean  normal  distribution  with  a specified 
standard  deviation  <rr,.  For  each  of  the  ranges  two  dis- 
tributions  of  the  measurement  errors  were  studied,  one 
with  standard  deviation  of  1 msec  and  another  with  10  msec. 
The  result  of  this  study  is  given  in  Table  1 and  the  reader 
should  refer  to  Fig.  El  for  parameter  definition. 

It  is  observed  that  both  range  and  bearing  estimates 
are  biased.  Since  the  output  of  the  position  estimator 
is  an  unbiased  estimate  there  are  doubts  as  to  the  random- 
ness of  the  time  perturbation  sequence  used  for  this  study. 
It  must  be  pointed  out  that  exactly  the  same  set  of  random 
numbers  was  used  for  each  range  and  measurement  error  case 
presented  in  the  Tables. 


The  ERIM  data  was  collected  at  a rate  of  300 

■ 

samples  per  second.  In  this  case  measurement  errors 
are  expected  to  be  at  least  of  the  order  of  3 msec. 
Inspection  of  Table  I reveals  that  under  these  cir- 
cumstances the  bearing  error  might  be  acceptable 
although  the  range  estimate  is  certainly  poor. 

In  order  to  obtain  a better  range  estimate  a 
larger  aperture  array  (than  the  ERIM  array)  is  needed 
at  the  ranges  considered.  Signal  reception  was  simulated 
to  an  array  consisting  of  18  geophones  formed  by  two  sub- 
arrays identical  to  the  ERIM  array,  but  separated  by 
1 km  (Fig.  E2) . Results  of  the  simulation  for  this 
array  are  given  by  Table  II.  The  improvement  in  range 
( estimation  with  this  array  is  significant  and  to  deploy 

such  an  array  in  the  field  should  not  add  considerably 
to  the  labor  involved  in  laying  a smaller  one.  An  in- 
vestigation on  optimal  array  configuration  is  planned. 
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Appendix  F 

Maximum  Entropy  Spectral  Estimation 


The  two  most  coition  methods  of  power  spectral  density 
estimation  are  the  periodogram  method,  in  which  the  spectrum 
is  estimated  by  either  averaging  over  different  samples  or 
smoothing  in  the  frequency  domain  the  Fourier  transform  of 
the  time-series  data,  and  the  Blackman-Tukey  approach.  In 
this  latter  method  the  autocorrelation  function  is  estimated 
from  the  data  using  lagged  products  with  a maximum  lag  taken 
typically  as  one  tenth  of  the  data  record.  The  Fourier 
transform  of  this  function,  after  proper  tapering  has  been 
applied,  yields  the  desired  spectral  estimates. 

These  conventional  methods  of  spectrum  estimation  pro- 
vide excellent  results  for  stationary  time  series  when  the 
record  length  is  large  compared  to  the  reciprocal  of  the 
lowest  frequency  of  interest,  but  they  have  certain  draw 
backs  when  it  conies  to  analyzing  non  stationary  or  transient 
phenomena,  such  as  the  signal  coupled  into  the  ground  by  the 
recoil  of  a gun. 

Both  the  periodogram  and  the  Blackman-Tukey  spectrum 
have  associated  with  them  window  functions  which  are  in- 
dependent of  the  data  or  the  properties  of  the  random 
process  which  is  analyzed.  The  window  function  relates 
the  average  estimated  spectrum  to  the  true  spectrum;  in 
fact  the  spectral  estimates  tend  to  be  the  convolution 
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in  the  frequency  domain  of  the  window  function  and  the 
true  spectrum.  Selection  of  a window  function  is  not 
a trivial  problem  since  good  resolution  and  statistical 
stability  must  be  traded  against  each  other  when  the 
amount  of  data  available  is  very  limited. 

A radically  different  nonlinear  method  has  been 
developed  to  estimate  spectra  with  increased  resolution. 
Interest  in  this  method  was  initiated  by  Burg  (cf.  J.P. 
Burg,  "Maximum  Entropy  Spectral  Analysis",  Ph.D.  disserta- 
tion, Stanford  University,  CA  1975)  who  used  the  term 
"maximum  entropy  method"  to  describe  an  algorithmic 
method  of  estimating  power  spectral  density  directly 
from  the  time  series  data  without  making  any  assumption 
about  the  characteristics  of  the  series  outside  the  ob- 
servation interval.  The  philosophy  of  the  maximum  entropy 
method  (MEM)  is  to  design  a unique  filter,  based  upon  the 
information  contained  within  the  available  data  record, 
such  that  when  applied  to  the  series  the  output  is  a white 
sequence.  Since  the  output  spectrum  is  a constant  it 
follows  that  the  spectrum  of  the  input  data  is  proportional 
to  the  inverse  of  the  transfer  function  of  the  filter. 
Entropy  of  a gaussian  process  is  defined  as 

H’jfc]  t*su\4r  o) 

~fi* 


where  fN  is  the  Nyquist  frequency  and  S(f)  the  spectral 
density.  For  a discrete  stationary  process  an  equivalent 
expression  for  entropy  is 
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where  C(N)  is  the  autocovariance  of  the  process  X,-  i 


C(N)* 


'f(o)  f (i)---.ffw) 
ffi)  fro>  - - f(N-j) 
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Assume  that  the  first  M+l  lags  of  the  autocovariance 
function  ^(0)j  f fi),  . . • , f(M)  are  exactly  known. 

The  idea  behind  the  MEM  is  to  determine  ffA'NJ),  . ••  ,f(N) 

in  such  a fashion  that  the  entropy  of  the  process  is 
maximized  at  each  step.  It  follows  then  that 
is  determined  by  maximizing  / C (M+l) j with  respect  to 
which  is  equivalent  to  setting 


ffi)  ffo) ffM-4) 
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Consider  now  that  the  process  can  be  described  as  an 
autoregressive  process  of  order  M; 


M) 


(5) 


where  €*n  is  a white  noise  series.  Multiplying  by  Xn„ 
and  taking  expected  values: 


f(*0  * «»  f (fc-O*  <*»  f • ■ ♦ o<M  f (fe-M) 


since  «0  f»r  *> o 

Substituting  k=l...,M+l  in  (6)  yields  a set  of 
equations  known  as  the  Yule-Walker  equations 

( fd)-  fto)  - - ‘X'm  f (M-ii  - O 
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Suppose  that  the  first  M+l  lags  of  the  covariance 
function  acre  known.  Prom  the  first  M equations  of  (7) 

the  coefficients  o/ljo<l) j c*H  may  be  determined 

and  the  unknown  ^ (M+l)  may  be  obtained  by  solving 

\ ?U)  f tO)  - - - - f(M-dj  I 
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Comparing  (4)  and  (8)  it  can  be  seen  that  MEM  spectral 
analysis  is  equivalent  to  fitting  an  autoregressive  model 
to  the  process. 

The  impulse  response  of  the  autoregressive  filter 
defined  in  (5)  is 

h (2)  - 1 - o^,  2 * £**■•***  “WmZ  ) 


where  z is  the  time  step  variable  and  the  power  spectrum 
of  the  linear  process  will  then  be: 


Appendix  G 

Adaptive  Autoregressive  Techniques 


It  was  proved  in  Appendix  F that  maximum  entropy 
spectral  estimation  may  be  obtained  through  fitting 
an  autoregressive  model  to  the  data  series; 

where  the  white  sequence  « 7 is 

termed  the  innovation.  In  the  MEM  developed  by  Burg 
the  filter  coefficients  are  determined  in  such  a fashion 
that  the  variance  Ge*  of  the  innovation  process  is 
minimized  over  the  observation  interval.  This  method 
provides  block  type  estimates  which  are  valid  over  the 
appropiate  design  segments. 

In  many  applications  where  not  only  the  existence 


but  the  evolution  in  frequency  of  signal  lines  is  to  be 
analyzed  an  adaptive  autoregressive  filter  is  in  order. 

A time  varying  whitening  filter  may  be  constructed  using 
a simple  algorithm  based  on  the  method  of  steepest 
descents.  The  filter  coefficients  are  updated  as  each 


data  sample  is  processed  in  the  following  fashion 
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and  &L  *s  a parameter  which  controls  the  adaptive  time 
constant  and  other  properties  of  the  algorithm.  Con- 
vergence of  the  algorithm  is  assured  provided  that 


M.* 


ra 

M <3? 


2 

where  is  the  input  power  level  and  0<fc<  2 
In  processing  the  ERIM  tapes  IAR  has  found  that  the 
algorithm  is  not  too  sensitive  to  variations  of  fa 
within  this  range.  The  adaptation  time  constant  of  the 
algorithm  expressed  in  number  of  data  points  is 


Finally,  the  instantaneous  frequency  spectrum  is 
given  by 
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technically  complete,  nevertheless  is  considered  very 
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